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We study a generalized double Jaynes-Cummings (JC) model where two entangled pairs of two- 
level atoms interact indirectly. We focus on the case where the cavities and the entangled pairs are 
uncorrelated. We show that there exist initial states of the qubit system so that two entangled pairs 
are available at all times. In particular, the minimum entanglement in the pairs as a function of the 
initial state is studied. Finally, we extend our findings to a model consisting of multi-mode atom- 
cavity interactions. We use a non-Markovian quantum state diffusion (QSD) equation to obtain 
the steady-state density matrix for the qubits. We show that the multi-mode model also displays 
dynamical preservation of entanglement. 
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I. INTRODUCTION 

The study of entanglement dynamics is crucial for the 
realization of quantum algorithms and quantum infor- 
mation processing protocols lj. A significant number 
of works have been devoted to study the dynamics of 
quantum entanglement under environmental effects 0- 
[7]. Many previous studies focused on the simplest sit- 
uation, namely, two-qubit entanglement dynamics. In 
reference [3], it was shown that contrary to what might 
be expected, the two-qubit entanglement can vanish com- 
pletely in a finite time which is often referred to as "en- 
tanglement sudden death" (ESD). One would naively ex- 
pect the two-qubit entanglement to decay asymptotically 
as a result of noise- induced decoherence effects [3] • ESD 
shows the fragility of entanglement under the unavoid- 
able interaction with the environment. A simple model 
for ESD is that of two two-level atoms interacting via 
Jaynes-Cummings (JC) Hamiltonians with two uncorre- 
lated single- mode cavities [S]. This double JC model is 
schematically depicted in Fig 1(a) Since the JC Hamil- 
tonians conserve the number of excitations (atomic plus 
photonic), the model can be treated as a four-qubit net- 
work. It turns out that even when both cavities are pre- 
pared in the vacuum state the entanglement between the 
atoms dies and revives periodically. This behavior may 
be interpreted as periodic entanglement transfer between 
atomic and photonic systems [5] . This model has also 
been extended to the multi-mode case where the atom- 
cavity couplings are described by a spectral distribution. 
For this model it was shown that entanglement cannot 
be protected regardless of the initial states (see [9l [10] 
and references therein). 

The purpose of this paper is to study the p reserv ation 
of entanglement in the network depicted in Fig |l(b) This 
model may be considered as an extension of the afore- 
mentioned double JC model. We assume that, initially, 



entanglement is only present in subsystems A1A2 and 
BiB 2 - In this network, subsystems Ai and Bi (i = 1,2) 
undergo excitation exchange interactions modeled via JC 
Hamiltonians. We restrict our attention to the situation 
where the cavities are prepared in the vacuum state. Un- 
der these assumptions, the model may be considered a 
four-qubit (atoms) and two-qutrit (cavities) system. In 
this context, we show that the pairs A1A2 and B\B 2 
can be prepared in certain partially entangled states such 
that they remain entangled at all times. This is the main 
result of this work. 

The paper is organized as follows. In Sec|H]we briefly 
discuss the double JC model and determine the corre- 
sponding evolution operator. In Sec jlllj we examine the 
entanglement dynamics in more complex scenario as por- 
trayed in Fig. (1(b)). We derived a compact expression 
for the evolution operator corresponding to the case of 
single-mode qubit-cavity interaction. This facilitates the 
study of entanglement dynamics for different initial states 
of the system. We show that there exist initial states for 
this network, such that the entanglement between the 
distant parties (A1A2 and B1B2) never vanishes. Addi- 
tionally, we studied numerically the minimum entangle- 
ment in these pairs as a function of their initial state. In 
this sense we found the optimal initial states of the pairs 
which, surprisingly, do not turn out to be maximally en- 
tangled. We also study the emergence of entanglement 
in the initially separable pairs A1B2 and A2B1. 

Finally, in Sec |IV| we include multimode atom-cavity 
interactions into our model. Here, the entanglement of 
the qubits is studied by means of a non-Markovian quan- 
tum state diffusion equations (QSD) [TlTfTS] . In addition, 
the residual entanglement in the qubits is determined in 
the steady state limit of the QSD equation The 
analytical results obtained in section corroborate the nu- 
merical results reported in Sec|III| 
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(a) Double JC model. Two two-level 
atoms, prepared in entangled state pa 
interact, locally, with uncorrelated single 
mode-cavities F\ and Fa. 




(b) Generalized double JC model. The 

qubits A\ and interact 
indirectly via common cavity modes. 

We assume that, initially, 
entanglement is only present in the 
pairs A1A2 and B\B2- 

FIG. 1. Double JC model (four-qubit model) and generalized 
double JC model (four-qubit-two-qutrit model). 



II. ENTANGLEMENT DYNAMICS IN THE 
DOUBLE JC MODEL 

Let the Hamiltonian acting on system (AiFi) be 



H, 



(i) 



H-nt, where 



(i) 



H ° 2 



and i = (1, 2). The spectrum of this Hamiltonian is well 
known [17] . The knowledge of the energy eigenstates and 
eigenvectors could be used to determine the time evolu- 
tion of the system. However, from a technical point of 
view, it is more convenient to find the time evolution op- 
erator by exponentiation of the Hamiltonian as de- 
scribed in |15j . It turns that this method may be also ap- 
plied to larger system as the one described in Fig. ( 1(b) I. 
For the double JC model we have: 



Ui := e 
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otNi 
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Here, we have assumed the zero detuning case (cu Ai = <^i) 



and used the relation [N 
show that 
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Clearly, the time evolution operator for the joint system 
AiA 2 FiF 2 is given U = Ui ® U 2 . Following [3|8], we as- 
sume that both cavities are initially in the vacuum state 
while the atoms start out in one of the following partially 
entangled states: 

\$a) = cos(a)\e Al ,e A2 ) +sm(a)\g Al ,g A2 ) (5) 
\V A ) = cos(a) \e Al ,g Aa ) + sin(a) \g Al ,e A2 ) . (6) 



Due to the fact that the JC Hamiltonian conserves the 
total number of excitations, the atomic reduced density 
matrix will be given by the X-state 



a / 

b e 

e* c 

\f* d 



(7) 



Throughout the present paper, we will quantify the 
entanglement E(p) by means of Wootter's concurrence 
C(p) [TB]. For states of the form Eq.|7|, the concurrence 
can be written in the compact form 



C(p) = 2max(0, |/| - Vbc, |e| - Vod). 



(8) 



Note that for the X-states of the form Eq.([5]) and Eq.Q 
we have C{p) = \ sin(2a)|. 

Using Eq.Q, one can determine the time-evolution of 
the reduced density matrix corresponding to the qubits 
A\A 2 8 . The entanglement dynamics is shown in 
Fig 2(a) and Fig |2(b)| Although both graphs describe 
the death and rebirth of entanglement, there are signif- 

In the 



(1) icant differences between Fig 2(a) and Fig 2(b) 



case where the atoms start out in the state \$ A ) given 
by Eq.|5]), the entanglement remains zero for finite pe- 
riods of time (except when a = 45°). These periods 
depend on the initial degree of entanglement in the sys- 
tem A1A2 (see Fig ]2(a) |. On the other hand, when the 
atoms are prepared in the state \^ A ) given by Eq.(|6|, 
the entanglement decays to zero periodically and recov- 
ers immediately, independently of a. Note that since we 
have assumed a symmetric scenario, the transformation 
a — > 7r/2 — a does not affect the concurrence. 



III. 



In 



A MODEL OF MUTUAL PRESERVATION 
OF ENTANGLEMENT 



this section, we study the scenario depicted in 
Fig. (1(b)). Here systems Ax, A%, B\, B2 are assumed to 
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(a) Concurrence as a function of time for system 
A1A2 when its initial state is given by |3>a)- The 
curves correspond to a = 60° (solid line), a = 45° 
(dashed line) and a = 30° (dotted line). 

E(At) 




(b) Concurrence as a function of time for system 
A1A2 when its initial state is given by • The 

curves correspond to a = 45° (solid line), a = 30° 
(dashed line) and a = 15° (dotted line). 

FIG. 2. Evolution of entanglement in the double JC model 



be two-level atoms while F\ and F2 represent single-mode 
cavities. Let the Hamiltonian acting on system (Aji^i^) 



be #W = H, 



H int where 



Following a similar method to that described in the 
previous section, we write the Hamiltonian ifW as 



H, 



I (A) " (B) h- t 

-oJAiCri " + -uj Bi al 11 + njjJiCii 1 a. 



H 



(i) 



and i=(l,2). The interaction of a single-mode quantized 
radiation held with N two-level atoms was hrst studied 
by Dicke [18 . The spectrum corresponding to the 
Hamiltonian Eqs. ([£])-( 10 1 was found long ago Q35] and 
its associated dynamics has been extensively studied in 
[20H22] . In addition, two- level atoms coupled to single- 
mode radiation held have been studied in connection 
with entanglement generation in cavity QED. (2"31 124) . 



ffW = huj t Ni + hX Ai c^ + h\ Bl C^ Bi \ 



with N t = ajai + Ua { z Ai) + a^>>) 
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CM = e B% a 



and e Bi 



(11) 
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From now on, we shall 



assume that the atoms and cavities are identical, that is, 
X Al = Xa 2 = As x = X B . 2 — X an d loi = 0J2 = oj. In 
addition, we shall again restrict our attention to the zero 
detuning case i.e., e A . = e B . = 0. These assumptions, 
plus the fact that [N^C^] = [N^C^] = 0, allow us 
to write the local evolution operator as 



iH (i) t 



—iutNi—iXtCi 
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where 



( <n a,i \ 

a\ a t 

a\ a, 

V a\ a\ J 



(14) 
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in the H Ai ®U Bi basis given by |lW) = |e-i, e^), (2^) = 
|ei,5i), |3 W ) = \g l7 e t ) and |4W) = \g ll g l ) . The operator 
Ui may be determined by exponentiating the matrix Cj. 
It can be shown that the even and odd powers of the 
operator Cj read: 



= 2" 
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where § 



Writing U; = 

e t^^-C^, we obtain the following compact 

expression for the evolution operator U,: 



4 
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(18) 



Clearly, the time evolution for the joint system 
A1A2B1B2F1F2 is given by Ui ® U2. We consider the 
situation where systems and B1B2 (systems A 

and B) are initially prepared in entangled pure states 
Pa = \4>a) (4>a\ and p B — \4>b) (4>b\- In addition, we as- 
sume that there are no additional correlations present in 
the total system. Thus, the initial density operator may 
be written as p = \4>a) (^a\ ® \<Pb) (4>b\ ® Pf x ® Pf 2 - At 
later times we have: 



p = Ui«)U2PoUl t «)U2 t . 



(19) 



Following the double JC model [7J [5] discussed in the 
previous section, we assume that the states \4>a) and \4>b) 
are of the form 



1$ 



MB) I 



cos(a) |e j4l (B 1 )' e ^2(s 2 )> 
sin(a)\9A 1 (B 1 ),9A 2 (B 2 )) (20) 



or 



I* 



MB)) 



cos(a) \e Al{Bl) ,g A2{B2) ) 
sin(a) |5Ai(fli), e A 2 (B 2 )) 



(21) 



In either case, we may write \4>a) = Sfe s fc l<?W,fe> <t>A 2 ,k)- 
Making use of equations Eq.( 18 ), Eq.(|19| and tracing out 
the degrees of freedom of systems BiB 2 and FiF 2 , we 
obtain the following expression for the reduced density 
matrix corresponding to qubits A\ and A2 (system A): 



Phi 



£^Tr B (p B V^ 



r (B 2 ) N 



(22) 



1,3 



Here 7r(l) = 1, 7r(2) = 2, for partially entangled states 
of the form Eq.Q while tt(1) = 2, vr(2) = 1, for par- 
tially entangled states of the form Eq.(|2"T|). The above 

V Bl(2) 
' ijkl 



operators are computed from the evolution oper- 



ator Eq.(18). They are given by 



^ k ? )] = ^f U2) {pf H2) m\ {2) \]) j 



(fc|Ul(2)|0 Al(2) ) 

(23) 



where |1) = |e) and |2) = \g). In the appendix Sec VI 



we list the set of non-vanishing operators VL fc j for the 
case where the cavities have a well defined number of 
excitations (i.e. pp i = \N) (N\). Note that expression 
Eq.([22| also holds true in the case where the systems 
A1A2 and B1B2 are prepared in different types of states, 
e.g., \$>a) ® \^b)- From symmetry considerations, we 



easily see that if the qubits start out in either |$a)®|^b) 
or \^ A ) ® then p A = p B at all times. 

Similarly, we can write down an expression for the re- 
duced density matrix for the qubit pairs AiB 2 and A 2 Bi. 
For simplicity, we shall consider only the situation when 
the initial state of the qubits is of the form \$a) <8> \®b) 
or \^ A ) ® Then we obtain: 



a 1 b 2 

r kl ran 



S * (*l V ir(j)nln(p) \<i) S q Sj (j\ V^(i) 



\p) V 



The usefulness of expressions Eq. ( 22 1 



(24) 

Eq.(23| and 



Eq.(24) lies in the fact that they can be evaluated in 
an automated fashion. They can also be applied to the 
more general case in which the cavities are prepared in 
mixed states [231. 



A. Partially entangled Bell States \$ A ) and |$s). 

We start by considering the case where systems A1A2 
and B1B2 are both initially in the same partially entan- 
gled state of the form Eq.(20|. As mentioned before, in 



the symmetric scenario in which the cavities are initially 
in the same quantum state, it suffices to compute the re- 
duced density matrix of one of the systems, say A±A2- 
Making use of Eq.(22| we determine the non- vanishing 
matrix elements 



A 2 4/ n b 2 + h 2 + 2p 2 . 2 . , 
pf x = a 2 cos 4 (a) H — sm 2 (2a) 

+ k 2 sin 4 (a) 

A , if \ , bf + hm-2p 2 . 2 
P22 = ad cos (at) H sm (2a) 

+ kn sin (a) 

A A 
P33 = P22 



p A = d 2 cos 4 (a) + 
+ n 2 sin (a) 



f 2 + 2p 2 + m 2 



sin 2 (2a) 



- A - " A - -e~™ {(<? + <?) ax? {a) 



Pu = Pa 



+ {I 2 + r 2 )sin 2 (a))sin(2a). 



(25) 



(26) 
(27) 



(28) 
(29) 



The functions a, b, c . . . can be found in the appendix 
Sec |VI| Of particular interest is the situation where 
both cavities are initially in the ground state, that is 
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It turns out that for certain values of a, the pairs A\A 2 
and B\B 2 remain entangled at all times (see Fig.( |3(a)] )). 
It is interesting to study the minimum entanglement 
E m in = min t C{p A {t)) in these pairs as a function of 
a. Based on numerical analysis, we conclude that for 
37.2° < a < 90°, there is always some residual entangle- 
ment in systems A\A 2 and B1B2, as shown in Fig.Q. 
This result is corroborated by Fig.( |3(b)"] ) where the time 
evolution of entanglement is shown for some values of 
a < 37°. If we adopt E min as a measure of the robust- 
ness of entanglement, we see from Fig.( |3(b)| ) that the 
most resilient state corresponds to a w 65.5° for which 
C(p A ) > 0.24. Interestingly, it does not correspond to 
the maximally entangled state (a — 45°). This non- 
trivial reflects a trade-off between the initial energy of 
the system and its entanglement. As a approaches 90°, 
the initial state of system approaches the energy eigen- 
state \g,g) <8> \g,g) <£> |0i) ® |0a). Consequently, its small 
amount of entanglement will not vary considerably with 
time. On the other hand, when a < 45°, the pairs (A1A2 
and B1B2) are more likely to be excited which renders 
dynamics of system more complex and tends to degrade 
the entanglement in the pairs. Moreover, the initial en- 
tanglement goes to zero as a approaches zero degrees. 
These two facts combined give rise to the critical value 
a cr 37 such that the entanglement in the pairs van- 
ishes for finite periods of time when a < a cr . 

It is also important to mention that in order to retain 
some entanglement in the pairs A\A 2 and BiB 2 , they 
must both be initially entangled. One can show that if 
the qubits start out in the state \$a) ® \q)b ® \q)b , 
then entanglement in A\A 2 will vanish for finite periods 
of time. Some entanglement will be transferred to B\B 2 
and for certain values of a it is possible to have one en- 
tangled pair at all times [25 . 

It is also interesting to look at the entanglement be- 
tween qubits A\ and B 2 (A 2 and B\). Here we shall also 
consider a symmetric configuration. Thus, it suffices to 
determine the density matrix for one of pairs, say A\B 2 . 
Making use of Eq. ( 24 ) we obtain 



p£ B » =a 2 cos 4 (a) + 

+ k 2 sin 4 (a) 
p£ B * = ad cos 4 (a) + 

+ fcnsin 4 (a) 

n MB 2 _ n A t B 2 
P33 — P22 



bh+p 2 . 2 



sin 2 (2a) 



fh + bm- 2p 2 



sin 2 (2a) 



p£ B * = d 2 cos 4 (a) + /T " 2 +p2 sin 2 (2a) 
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+ lr sin 2 (a)) sin(2a), 
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(34) 
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(a) Concurrence as a function of time for the cases 
a = 45° (solid line), a = 60° (dashed line) and 
a = 85° (dotted line). 
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(b) Concurrence as a function of time for the cases 
a = 37° (solid line) and a = 15° (dashed line). 

FIG. 3. Evolution of entanglement for systems A1A2 (B1B2) 
The qubits are initially prepared in CS> |$s) ■ 
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FIG. 4. Minimum Entanglement in A1A2 (B1B2) as a func- 
tion of Q. 



where the functions a, 6, c . . . can be found in the ap- 
pendix Sec|VI| 

Note that the pairs A^B 2 and A 2 B^ start out in separa- 
ble (mixed) states. As a result of indirect interactions be- 
tween AiBi, some fraction of the original entanglement in 
AiA 2 and B\B 2 will be transferred to these pairs. We as- 
sume that the cavities are prepared in the vacuum state. 
The case where A\A 2 and B\B 2 are initially in the state 
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$ B ) with a — 75° is shown in Fig.(|5(a)[). From B. Partially entangled Bell States and \^ B ) 



this graph we see that for certain periods of time, we 
have the choice of selecting either two entangled or two 
separable pairs. Note that for this value of a, the con- 
currences exhibit an approximately sinusoidal behavior. 
The dynamics corresponding to a = 45° turns out to be 
far more complex as shown in Fig.( |5(b)"j ). Note that these 
graphs suggest that the entanglement in A\B<i can never 
exceed that in B1B2 (or A1A2) ■ In fact, this is a direct 



consequence of equations Eqs.(26), (291, (31) and ( 34 1 



Using the inequalities cq < \{c 2 + q 2 ) and It < 1 ^ r one 
proves that Ipf^ 52 ! < \Pu\- In addition, we have 



It turns out that if systems A1A2 and B1B2 are ini- 
tially in states of the form in Eq.(21), their entanglement 



cannot be preserved, for any value of a. Moreover, as in 
the previous subsection, the concurrence of A1B2 (A2B1) 
never exceeds that of A1A2 {BiB2). Using expression 
Eqs.(22| and (23), we compute the density matrix de- 



scribing A1A2 and B\B2- The nonvanishing matrix ele- 
ments read 



P22 P22 2 



(h-b)(m-f) 2 
sin (2a) 



~cos 2 (2\t\l N + sin 2 (2a) < 



(35) 



which completes the proof. 
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(a) Concurrence as a function of time for A1B2 
(dashed line) and A1A2 (solid lines). Here a = 75 
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(b) Concurrence as a function of time for A1B2 
(dashed line) and A1A2 (solid lines). Here a = 45°. 
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Pi X = afccos 4 (a) H sin 2 (2a) + ak sin 4 (a)(36) 



P22 = an cos 4 (a) + 

+ dfcsin 4 (a;) 
p^ 3 = dk cos 4 (a) + 
an sin 4 (a) 



fh + bm- 2p 2 
4 



sin 2 (2a) 



A J7 4/i , fh + bm-2p 2 . 2 

" ' — sin (2a) 



4 

fm+p 2 . 2 



P4 4 = dn cos 4 (a) + - "" ' - sin 2 (2a) 

+ dn sin 4 (a) 

A A * cl + qr 
P23 = P32 = — o — sm(2a). 
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The curves corresponding to entanglement as a func- 
tion of time for different values of a are shown in 



Fig. (6(a)). We found that the concurrence vanishes for 
finite periods of time for every a. Note that the curves 
appear to coalesce and go to zero simultaneously regard- 
less of the value a. However if we zoom in on a portion 
of the graph, we can see that this is not the case (see 
Fig.( |6(b)"| )). We conclude this section with the following 
remark. Expressions Eqs.(22) and (23) may also be ap- 



plied to the situation when the cavity modes are excited, 
i.e. N > 1. We found that entanglement cannot be pre- 
served (in the sense of Fig 3(a) ) unless both cavities are 
prepared in the vacuum state [N = 0). The same holds 
true for the \^a) ® \^b) case. 



IV. MULTIMODE INTERACTION 



FIG. 5. Evolution of entanglement for systems A1A2 (B1B2) 
(solid line) and A1B2 ( A2B1) (dashed line). The qubits are 
initially prepared in \&a) <S> |$b) ■ 



In this section we extend our analysis to the case where 
the two-level atoms interact with multi-mode cavities 
(structured environment). We describe this situation 
with the following natural generalization of the Hamil- 
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reduced density matrix is constructed from 




10 12 14 



(a) Concurrence as a function of time for the cases 
a = 25° (solid line), a = 45° (dashed line) and 
a = 60° (dotted line). 



p = m[\m*'))(M**)\] = 



IT 

(46) 

Since there is no direct interaction between the sub- 
systems (Ai,Bi,Fi) and (A 2 ,B 2 ,F 2 ), the noises gener- 
ated by the two local baths and O-operators are uncor- 
rected. The zero temperature assumption together with 
the chain rule = Jl* ds^if- j^r, allow us to construct 
the following QSD equation 



E(At) 
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(b) Concurrence as a function of time (zoomed in) for 
the cases a = 25° (solid line) and a = 45° (dashed 
line). 

FIG. 6. Evolution of entanglement for systems A\Ai (B1B2) 
The qubits are initially prepared in \^a) <8> |*s) • 



tonian Eq.Q and Eq.([To|): 



rrO) _ Tj(i) 1 rr(i) 
H - AB + H F 



H 



(i) 
AB 

r(0 



hujj 



{At) 



r(Bi) 



Hp = ^2 hwika-ik' a,ik 
fc 



(41) 
(42) 
(43) 

(44) 



Here, aife and a ik , correspond to the creation and annihi- 
lation operators of the k th electromagnetic mode in the 
i th cavity and frequency Wik- Using the Bargmann state 
of the baths to trace out the baths' degree of freedom, 
we have 



id t Mz*)=Y,( H AB+^l>2 X *K 



fee 



-iuiikt 



dz 



(45) 



ik 



where ipt(z*) is the system stochastic vector for the four 
qubit system A\A 2 B\B 2 and Li = 



{At) 



The 



L\ I dsGj(t, s)Oi(t,s,z*)}i> t (z*), 



n- 



^AB + LiZ* t 



(47) 



where z* t = -i£ A A* fe z* fc e-« \ G&a) 

tt |2 e -*-.*(*-) and O t (t,s,z*)A(z*) = s^M**) 



F il (t)0 1 +F i2 (t)0 2 - 



The QSD method yields O l (t, z*) 
iUi(t, z*)0 3 , where 0\ = L, 2 
3 = (t (A) (t (B \ and Uifaz?) = f*dsUi(t,s)zl to be 



ai A) a iB) 



*i A) oi B \ 



determined. Assuming Gj(i, s) 



Eli 



-H\ts\ which 



corresponds to the Lorentz spectral density for the 

1 r 7| 



multi-mode cavities S(u)j 

r 7; 



, we obtain 



dtFn(t) 



+ (- 7i + tui)F a + + 3F, 



1 2 



-2 U » 

d t F l2 (t) = (-7j + Wi)F a - Fl + 4F a F a + F% 
-ft' 

d t Ui(t) = -2i 7l F t2 + (-2 7i + 2tUi)Ui 



(48) 



(49) 



(50) 



where Ui(t) = J Q dsGi(t, s)Ui(t, s). The boundary con- 
ditions are given by -Fii(O) = F i2 (0) — f/,-(0) = and 
Ui(t,t) = -UF Q {t). 

It is also known that an open system in a non- 
Markovian bath can approach a stable final state in the 
long time limit, as long as the bath correlation function 
has a well-defined Markov limit. For our model we have 
that 7i — > 00 implies G(t,s) — > TS(t,s). For notational 
simplicity, we write the density matrix for the total sys- 
tem as 



(51) 




where a,b,d,- ■ ■ ,p represent a 4 x 4 sub-matrices. We 
work in the basis {|e, e, e, e) , \e,g, e, e) . . . \g, g, g,g)} € 
Ha! ®'Hb 1 ® Ha 2 ® %b 2 ■ The final stable state for this 
model is found to be 



/ 

F G H 

J K L 

V N O P 



where 



F = K=-G = -J, H = -L = N^ = -0 ] , 
and each non-vanish sub-matrix has the form 



/ 

v —v 

— v v 

V w* -io* 



(52) 



(53) 



(54) 



Just as in the previous section, we shall consider the case 
where the four qubits are initially prepared in a state of 
the form ® |$ B ) or \V A ) ® (see Eqs.Qn}). 

For these two cases, tracing out the degrees of freedom 
corresponding to any pair of qubits we obtain the reduced 
density matrix for the other pair. Thus, in the long time 
limit we find that 



(55) 



and 



,,AiB 2 _ 



F22 


Q 





H24 





F33 














K22 





N 42 








K33 + P44 


F33 








H34 





F22 














#33 





N 43 








K22 + P44 



(56) 



In particular, for the case where the qubits are initially 
in \&a) ® \&b), we obtain 



y 








a; 





y 














y 





X* 








1 - 



(57) 



3y 



and 



A 1 B 2 _ 



1 


l> 





-x 







/; 













v 


\ 


-as* 





1-3?/ 



(58) 



where y = 1/4 cos 2 (a) sin 2 (a) and |x| 
l/2cos(a) sin 3 (a). From the above expressions of 



Eq. (|57J) and Eq. (58 1, we conclude that in the long time 
limit we have C(p^ B2 ) = C(p AlA2 ). The concurrence 
is given by C — 2max{(|a;| — y),0}. Note that the con- 
currence does not vanish for partially entangled states 



having arctan(0.5) w 26.6° < a ^ 90°. Interestingly, the 
maximum of the concurrence C max = 0.24 is attained 
at a R3 65.3° which is consistent with the single mode 
model discussed in section Sec lIIII 

As for the case ( \^a) ® |^b)), w e find that the long 
time density matrix now reads 



/ y 

y 

y 

V l-3y 



(59) 



Therefore there is no entanglement present in the final 



state of A1A2. Following the similar steps, p 



AiB 3 



also 



ends up as a separable state for this initial condition. 



V. CONCLUSIONS 

In this paper we studied the entanglement dynam- 
ics in a generalized double JC model. We showed that 
although the system evolves non-trivially, two pairs of 
qubits (A1A2 and B1B2) can preserve some fraction of 
their initial entanglement. We found a family of initial 
states for the system \$ A ) ® (8) |0i) ® |0 2 ). such 
that the entanglement in the pair never vanishes. We 
also determined the optimal initial state for which the 
concurrence is greater than 0.24 at all times. Interest- 
ingly, this optimal state is not a maximally entangled 
state. This result does not involve conditional dynamics 
(i.e. no quantum measurements are required). The sce- 
nario presented in this paper should be compared with 
the double JC model (see Secjll]) where this preservation 
is not possible for any initial configuration of the system. 
Thus, putting aside questions related to the experimental 
realization of our scenario, the comparison of both mod- 
els suggest that storing the qubits in pairs may be a way 
to protect their entanglement. One can envision even 
larger networks with qubits prepared in multi-particle 
entangled states. It would be interesting to explore such 
systems and study the amount of entanglement available 
at any time. The aforementioned effect of mutual preser- 
vation can be interpreted as the result of the constructive 
interference of the amplitudes corresponding to processes 
of emission, absorption etc. It may be also interpreted 
as partial entanglement transfer, that is, the initial en- 
tanglement cannot be completely redistributed over the 
rest of the pairs. In the double JC model the initial en- 
tanglement of the pairs can be completely transferred to 
the cavities [5]. 

Naturally, one is tempted to study all pairwise quan- 
tum correlation and attempt to establish entanglement 
conservation rules for this model (as in [55] and |27j). 
For mixed states we only know the separability criteria 
for the low dimensional Hilbert spaces C M x C N with 
M = 2 and N = 2 or N = 3 )28| or for the case of 
bipartite Gaussian states j^Hj- As a result, all pairwise 



concurrences can be computed except for F1-F2 (cavity- 
cavity) which is, effectively, a 3 x 3 system. 
Finally in Scc |IV| we included multimode qubit-cavity in- 
teractions and studied the dynamics of the system by 
means of a non-Markovian state diffusion equation. We 
found the density matrices in the long time limit. The 
results corroborate those from the single-mode interac- 
tion model. 

The latter suggests that it would be interesting to 
explore other multi-qubit configurations and interaction 
models. Such studies may lead to a better understanding 
of the entanglement dynamics and provide interesting in- 
sights into the problem of protecting entanglement from 
the environment. 
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VI. APPENDIX 



The operators Yijki operators defined in Eq.( 23 ) satisfy 
the properties V\ jkl = Vi kj i and Y^k^ikkl = When 
the cavities are prepared in the pure state = \N) (N\ 
these assume the form 



Vim 


= (N 


V U 21 


= {N 


V1221 


= {N 


V2112 


= {N 


V2122 


= (N 


V2222 


= (N 


V1112 


= (N 


V1122 


= (N 


V1222 


= {N 


where 7 


— Ldt. 






V2212 = ^ 


' 2122' 



(elU^le) (e|Uj|e) \N) ■ 
{e\Vj\e) (g\Ui\e)\N) : 
(e\Uf\g) (fl|Ui|e)|JV) 
(g\Vj\e) (e\Ui\9)\N) 
(g\Uj\e) (<?|I%)|A) 
(g\uS\g) (g\Ui\g)\N) 
(e\Uf\e) (e\Ui\g)\N) 
(e\vS\e) (g\V t \g)\N) 
(e\Uf\g) (g\V t \g)\N) 
Additionally we ha\ 

^2211 = V f 1122 , V2111 



a 
b 



(60) 
qW(61) 
(62) 
(63) 
e* 7 (64) 
(65) 
(66) 
(67) 
(68) 



c 


d 

/ 

h 

k 

I 



m 

n 



p 

q 

r 





-p 

1211 = 

4ll2' V222I 



1121' 



^1222 an d V1212 = V2121 = 0, which completes the list. 



The functions a,b,c. . , read 



(1- 



N + l 
A + 3/2 



sin 2 (A<V^ + 3/2)) 2 



N + 1 



4(A + 3/2) 

„4 



sin 



'(2Xty/N + 3/2) 



b = cos i {Xt y / N + 1/2) 



A 



4(A+ 1/2) 
N+l 



sm 2 (2Xty/N + 1/2) 



(69) 



(70) 



4^(^+1)2-1/4 



sm 



(2Xt y / N + 1/2) sin(2Xty/N + 3/2) 



sin 2 (Xty/N + 1/2)(1 



N + l 
A + 3/2 



3in 2 (AV^ + 3/2)) (71) 



N + l 



(2Xty/N + 3/2) 



f = 
h = 
k = 
+ 
I = 



4(A + 3/2) 

sin 4 (AiV^ + l/2) + ,, A jV+ 1 1 /0 , sin 2 (2At^ + l/2)(73) 



sm 



4 (Xt^N + 1/2) 



4(A + l/2) 
A 

4(A + l/2) 



sm' 



: (2Atv'^' + V 2 )(74) 



A 



4(A- 1/2) 



sin 2 (2AtV^- I/ 2 ) 



^^sin^(AtV^l72) 
(A - 1/2) 2 V V 1 ' 



(75) 



A 



4(y/N* - 1/4) 



sm 



{2\tyjN - 1/2) sin(2A<v/A + 1/2) 



- sin 2 (Ai v /]vTl72)(l - 
to = cos 4 (A<v/A + 1/2) 



N 

N - 1/2 



sm' 



, JV + 1 , x sin 2 (2A<v/A + l/2) 
4(A+l/2) 1 V -r / ; 

A 

4(A- 1/2) 

a- w 



'(Xty/N -1/2)) (76) 



(77) 



sm' 



: (2A<v/A- 1/2) 



P : 
1 : 



A - 1/2 
1 

8(A + l/2) 



sm 



sm' 



\Xt^N -1/2)) 2 
! (2AV^ + l/2) 



(78) 
(79) 



cos 2 (Atv^vTT72)(l - sin 2 (A< v ^vT372)) 



A + 1 
4 v /(A + l)2-l/4 t 
sin(2Atv/A + 3/2) 



1 (2Xty/N + l/2) 



(80) 



r = cos' 



; (Atv/A + l/2)(l- 



A 

A - 1/2 



sm 



: (A<v/^-l/2)) 



A 



- 1/4 



sm 



{2XtyjN - 1/2) sin(2At V / ^ + 1/2) (81) 
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